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Abstract. In this study, wave equations with initial boundary conditions have been studied. The general form of 
the wave equation has been derived. The first order and second order difference schemes were established for the 
presented IBVP. The stability of the difference schemes has been guaranteed. The approximation solution of the 
problem was achieved by using finite difference methods. Two different examples are provided. A comparison 
between the exact and approximation solution has been carried out. Absolute errors of the problem have been 
presented by using MATLAB software. Moreover, the comparison shows that the second order difference scheme 
is amore accurate result than the first order. It is shown that the results of the comparison guaranty the reliability 
and accuracy of the presented method. 
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1. INTRODUCTION 


In the last decades, Partial differential equations (PDE) have focused on many stud- 

ies because of their common appearance in several applications in mathematics, physics, 
seismology, science, finance, engineering, and mechanics [1, 2]. PDEs are important for 
modeling a wide variety of problems in science and engineering, including electrodynam- 
ics, elasticity, wave propagation, signal analysis, and thermodynamics [3-5]. 
However, the exact solution of some types of PDE is known [1-3]. Much attention has 
been given to studies to test the reliability and accuracy of the approximation and numeri- 
cal techniques such as finite difference methods [4—6]. The wave equation is an important 
second order linear PDE for describing waves, such as mechanical waves or light waves 
[1-5]. It 1s the simplest example of the hyperbolic differential equation. Wave equations 
and modifications play fundamental roles in continuum mechanics, quantum mechanics, 
plasma physics, general relativity, geophysics, and many other scientific and technical dis- 
ciplines [3-5]. Nytrebych et al. presented an analytical method for studying the wave 
process mathematical model under two-point time conditions[6]. 
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Moreover, several different approximation methods were presented for solving wave equa- 
tions [7]. Finite difference approximations have been used for solving different types of 
partial differential equations in the works [8—15]. Qian and Weiss worked on Wavelets and 
the numerical solution of partial differential equations [11]. Ashyralev and Modanli pre- 
sented a finite difference method for finding a numerical solution for telegraph PDE [17]. 
Sweilam, et al. used the finite difference method to find numerical solutions of two-sided 
space-fractional wave equations [18]. Carcione and Helle presented a Numerical solution 
of the wave equation on a staggered mesh [19]. Saadatmandi and Dehghan worked on 
a onedimensional wave equation with an integral condition [20, 21]. On the other hand, 
various studies on stability estimates have been done for different types of PDEs [21-25]. 
However, different numerical methods are presented for solving wave equations with ini- 
tial boundary conditions [19-22]. In this work, we will present a numerical method based 
on difference scheme method for solving wave equations. 


2. DERIVING THE 1D WAVE EQUATION 


In this section, we will present the derivation of wave equations which have been done 
in [25]. The derivation of the one-dimensional wave equation starts from Newtons and 
Hookes law. As shown in Figure 1. the key notion is that the restoring force due to tension 
on the string will be proportional to the curvature at the point. Then the force, ku,,. should 
be equal to mass times acceleration Ppu;;, where c = \/k/p turns out to be the velocity of 
propagation [25], we obtain 


a: 
Uty = C Uxy 





FIGURE |. The restoring forces on a vibrating string, proportional to curvature. 


Let u = u(x,t) denote the displacement of a string from the neutral position u = 0. Let 
k = k(x), be the elasticity, and p = p(x) be the mass density of the string. Especially, We 
do not assume the string is uniform in this derivation. Consider a small piece of string in 
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the interval |x, x-+ Ax]. Its velocity be u;(x,t), itis mass p(x)Ax, and thus its kinetic energy 
iS: 


The total kinetic energy of the string is provided by an integral, 


1 sel ) 
K=5 p.(uz) dx. 
2 Jo 


The potential energy for a string is (5)y given from Hookes law, where y is the length of 


the spring [25]. The length of the string is presented by arclength, in the case of a stretched 
string, as follow 


ds = \/1+uzdx, 


and potential energy is given in the following form 


Lk 5 
cy ~(1+u%)dx. 
0 2 


From it, the function u is presented as the integral over time of the difference of these two 
energies, SO 


fa SL [pay -elt-+ uP ae 


Adding 6 times a perturbation 4 = h(x,t) to the function u, we obtain 
T pl 
L(u+ 6h) = L(u) + 5 | / p.u;.h, — k.u,h,dxdt + highr order in 6. 
0 JO 


The principle of least action says that for u to be a physical solution, the first-order term 
should vanish for any perturbation h. By applying integration by parts, we obtain 


/iese 
0 -| / (—P Upp + KU + ky) dxdt. 
0 JO 


Because of all choices of h, this integral 1s zero, the first factor in the integral must be zero, 
then we get the wave equation for an inhomogeneous medium, 


P Ug = Ky, + ky.tty 


above equation reduces to the usual two-term wave equation, when the elasticity k is con- 
stant, 


a 2 
Utt — C Uxx 


[k 
where the velocity c = a varies for changing density [25]. 
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By applying initial boundary conditions, we can state the problem as follow: 
gt C6 Ut *) On x1, 
u(t,0) = u(t, z) = 0, Of <T, (1) 
u(0,x) = w(x), u,(0,x) = (x). 


Where w(x), P(x) (x € [0, 2]) and f(t,x) =0, ((t,x) € [0,1] x [0, z]) are smooth functions, 
problem (1) presents a wave equation [11]. 


3. STABILITY AND DIFFERENCE SCHEME 


This section introduces the explicit and implicit methods for approximate wave equation 
solutions in (1). Moreover, von Neumann’s Stability Analysis will be presented. A Simple 
Explicit Method can be applied to (1), by using central differences for both time and space 
derivatives, we obtain 


Z k kak 
Unt | — 2, + un | i a 


n—1 
2 yp (2) 
CT ae 
If we set a= 7° we can simplify (2), we get 


uk) — —yk-1 491 —a)uk 4.02 (uh +uh_;] (3) 


The stability of the explicit scheme in (3) guaranteed in [24—26], we state it with the usual 
ansatz [25], 
ek — gkgiktn 
From which we obtain the following expression for the amplification factor g(k) 
g’ = 2(1—a’)g—1+4+2a’g cos(kAx). 


After several transformations, the last expression becomes just a quadratic equation for g, 
namely 


g —2Bg+1=0, (4) 
where B = 1 — 2a’sin? (>) . By using quadratic rule solutions of the equation for g(k) 
1S 
g12=B+VB-1. 


Notice that if |B] > 1 then at least one of the absolute values of g;2 is bigger than one. 
Therefore one should desire for B| < 1, we get: 


g12=B Fiv Be —1. 
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and |g|? = B* + 1— B? = 1. 
That is, scheme (3) is conditional stable. The stability condition read 


kh 
-1< | —2a*sin*(>) <1 


; ; - cT 
what is equivalent to the standard CFL condition a = — < 1. 


For the wave equation in (1), we will introduce an implicit scheme from (2). However, 
the second-order difference scheme is used for presenting the conditional stability of the 
problem. By replacing all terms on the right-hand side of (2) by an average from the values 
to the time steps k — 1 and k+ 1, we obtain 
Wht Duk uk! 
72 ~ D2 (uk 
The difference equation in (5) is called a second order difference scheme for the problem 
(1). We used the implicit method. The stability of the difference scheme (5) is guaranteed 
in [25-28]. We will state it is stability, by using the standard ansatz 


— py ae ty + uw — 2uk—! + uh) (5) 


n—| 


gk — gk giktn, 
Leading to the equation for g(k) 
Bg —2¢+B =0, 
where B = 1 +2a?sin? (F). We have B > 1 for all k. Hence the solutions g; 2 take the 


1+i,/1—- Bp? 


£1,2 = : 
p 


form 


and 


_ (1 — R2 


That is, the implicit scheme (5) is absolute stable. 


4. THE METHOD 


Both difference equations in (2) and (5) can be rewritten to the following matrices 
formulas are obtained as 
AU 41 + BU, + CU,_1 = Don, l<n<M-1, u=uy=0. (6) 


where, A,B and C are (V+ 1) x (N+ 1) matrix, Uj+1,Un,Un—1 and @, is (N+1) x 1 
vectors, in the present work @, = 0 [11-13]. We applied a modified Gauss elimination 
method to solve the difference equation in (6). Then, we are looking for a matrix equation 
solution, which was presented as follows [16]: 


Uj = Oj41Uj41 + Pj41; uy =0; j=M—1,...,2,1. 
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Where 8; are (V+ 1) x 1 column vectors, and @; are (V+ 1) x (N+ 1) square matrices, 
defined by 
j41 = —(B+Ca;) ‘A, 
and 
Bj.1=(B+Ca;)'(De—CB,), j=1,2,...,M—1, 

where j = 1,2,...,M—1, By, is the (V+1) x 1 zero column vector, and a, is the (N + 
1) x (N+ 1) zero matrix. The Matlab program computes the results. For three distinct 
values of M and N, the comparison was made between numerical and exact solutions. The 


; _— 7 1 
maximum error indicated where h = WT and T = N’ Maximum absolute error calculated 
by 


EN = max |u(tg,x,)) —ué 


1<k<N-1 
1<n<M—-1 








where exact solution denoted by u(t,,x,) and approximation solution denoted by uk at 
points (t,,Xn). 


5. NUMERICAL COMPUTATION 


This section will apply the first order difference scheme in (2) and the second order 
difference scheme in (5) for two different numerical examples for the wave PDE with 
initial-boundary conditions. For both examples, obtained exact solutions will be compared 
with the numerical result for testing the method’s ability. In applications, let us consider 
the general form of the IBVP hyperbolic wave partial differential equation 


Ur (t,X) = C7Uxx(t,x) O<x<L, 
u(t,0) = u(t,z) =0, Os taT, 
u(0,x) = W(x), u(0,x) = (x). 

Example 5.1. Consider the following IBVP for the wave equation 


Uy (t,X) = Uxx (t,x) O<x<lL, 
u(t,0) = u(t, az) = 0, O<t, (7) 
u(0,x) = sin(x), u;(0,x) =0. 


Analytical solution of the problem (7) can be found by using one of the exact methods, 
which is u(t,x) = cos(t)sin(x). Then the first order difference scheme for (7) is presented 


as follow: : ; : 
he _ 2ui + 7 i __ Unt) — QU; “ap Uy 


7 7 h 
uj =u,=0, O<n<NO<m<M (8) 

0 
u,—u 
" 0 
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For the second order difference scheme, we obtain 





k+1 k k—1 k+1 _ 94,,k+1 k+1 k-1 _ 4,,k-1 k—-1 
Uy — 2u, + uy, = Uns} 2Uy, see TU y 2Uy, TU, 4 
72 h2 
us=uk=0, O<n<N,O<m<M (9) 
1 0 2 1 0 
un —Uu Tue —2u.+u 
O_ @: n n __ n n n 
LSI iy) = - 4 * 


T 2 T 


The error analysis of the problem (7) is shown in table 1. However, the figures of exact 
and approximation solutions have been shown in figure 2. 


Example 5.2. Investigate the following IBVP for wave partial differential equation 


lt 2) g(t) O0<x<L, 
u(t,0) =u(t,) =0, 0<r, (10) 
u(0,x) = 4sin(x) — sin(2x) — 3sin(5x), u,(0,x) = 0. 


The exact solution of the problem is u(x,t) = 4cos(t)sin(x) —cos(2t)sin(2x) —3cos(5t)sin(5x). 
The first order difference scheme for the problem (10) is 


Unt! = 2un tun! Ung ~ 2H Fn 

a Le ee 

i =1, =, 0<n<N,0<m<M (11) 
1 0 

u2 = 4sin(x,) — sin(2x,) — 3sin(5xp) ac sl 

nn n n n/» t Ts 


The second order difference scheme (Implicit method) can be obtained as follow: 


k+1 k k-1 k+1 k+1 k+1 k-1 k-1 k-1 
Uy 2p Fn Mn An En tM nF Mn 


n—1 
7 7 h? 
uj=uk=0, O<n<N,0<m<M (12) 
1_ 10 3 1,0 
= Ci = 2 
u2 = 4sin(Xp) — sin(2x,) —3sin(5Xxn), es ee eee 


2 T 


The error analyzes of problem (10) is shown in table 2. However, the figures of exact and 
approximation solutions have been shown in figure 3. 
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TABLE 1. Error analysis for exact and approximation solution of Ex- 
ample 5.1 


M N 


Difference scheme in (8) | 1.727910~2 | 8.53731072 | 4.23821073 


Difference scheme in (9) | 8.448210~* | 2.1015107* | 5.234110-> 





TABLE 2. Error analysis for exact and approximation solution of Ex- 
ample 5.2 


Difference scheme in (11) | 5.5646107! | 2.3548107! | 1.0444107! 


Difference scheme in (12) | 3.4757107! | 9.122210~2 | 2.277110~2 





The exact and approximation solutions of given examples are presented in the following 
figures. 


Exact solution Approximate solution 











FIGURE 2. Exact solution and approximation solution of the problem 
(7), where N=M=25. 
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Exact solution Approximate solution 
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FIGURE 3. Exact solution and approximation solution of the problem 
(10), where N=M=50. 


6. CONCLUSION 


In the presented article, the wave equation was modeled mathematically, and the prob- 
lem was verified. The initial boundary value problem for the wave equation has been 
discussed. The difference schemes of the first order and second order of accuracy for wave 
partial differential equations are showed. Stability for the problem and difference scheme 
has been established. Modified Gauss elimination technique based on difference equa- 
tions used for achieving approximation results. An algorithm has been established, and 
the programs are written by MATLAB software for numerical computation. Two different 
examples with the analytical solution are given. Approximation results compared with the 
exact solution. Good results were obtained for the problem. Error analysis tables have 
been presented, and it is guaranteed the accuracy of the method. The comparison shows 
that the second order difference scheme is a more accurate result than the first order. It 
is worth mentioning that the technique can be used as a very accurate algorithm for the 
presented type of wave PDE. The accuracy of the numerical results supports theoretical 
statements of the presented solution of finite difference methods. 
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